01_airy

Uvod¶

Airyjevi funkciji $\text{Ai}$ in $\text{Bi}$ se v fiziki pojavljata predvsem v optiki in kvantni mehaniki. Definirani sta kot neodvisni rešitvi enačbe

$$ y''(x) -xy(x) = 0, $$

ki je znana kot Airyjeva oz. Stokesova enačba. To je najenostavnejša linearna diferencialna enačba drugega reda z obratno točko (tj. točko, kjer se značaj rešitev spremeni iz oscilatornega v eksponentni). Neodvisni rešitvi enačbe sta predstavljivi v integralski obliki kot

$$ \text{Ai}(x) = \frac{1}{\pi} \int_0^\infty \cos (t^3/3 + x t) \,\mathrm{d} t \>,\quad $$

$$ \text{Bi}(x) = \frac{1}{\pi} \int_0^\infty \left[ \mathrm{e}^{-t^3/3 + x t} + \sin (t^3/3 + x t) \right] \,\mathrm{d} t \>. $$

Naloga¶

Z uporabo kombinacije Maclaurinove vrste in asimptotskega razvoja poišči čim učinkovitejši postopek za izračun vrednosti Airyjevih funkcij $\text{Ai}$ in $\text{Bi}$ na vsej realni osi z absolutno napako, manjšo od $10^{-10}$. Enako naredi tudi z relativno napako in ugotovi, ali je tudi pri le-tej dosegljiva natančnost, manjša od $10^{-10}$. Pri oceni napak si pomagaj s programi, ki znajo računati s poljubno natančnostjo.

Lastnosti Airyjevih funkcij¶

Funkciji $\text{Ai}$ in $\text{Bi}$ sta definirani na celotni realni osi. Za vrednosti $ x \in \mathbb{R}^-$ sta obe funkciji oscilatorni in omejeni. Za $ x \in \mathbb{R}^+ \cup \{0\}$ pa je funkcija $\text{Ai}$ omejena, medtem ko $\text{Bi}$ narašča čez vse meje.

Airy Functions

Graf funkcij $\text{Ai}$ in $\text{Bi}$

Matematični pristop in potencialne težave¶

Za izračun vrednosti funkcij integralov ne bomo neposredno izvrednotili, temveč bomo funkcije razvili v Maclaurinovo in asimptotsko vrsto. Pri tem bomo pazili, da formule implementiramo rekurzivno in tako čim bolj zmanjšamo število potrebnih računskih operacij. Na grobo ocenimo, da z uporabo rekurzivnih formul zmanjšamo časovno zahtevnost programa iz $O(n^2)$ na $O(n)$, kar v računalniškem lingu razume kot "iz zelo počasi v precej hitro".

Za natančnost bomo poskrbeli z ustreznim izborom števila členov v Maclaurinovi oz. asimptotski vrsti. Ker ena vrsta podaja dober približek za zelo majhne vrednosti $x$, druga pa za zelo velike, bomo vrste v točki $x_0$ zlepili. Točko $x_0$ bomo določili na podlagi števila členov, ki jih potrebuja ena in druga vrsta za doseganje tolerirane napake.

Potencialne težave bodo (in tudi so) predstavljale predvsem: omejitve podatkovnega tipa float, iskanje zanesljive referenčne vrednosti za oceno napak in iskanje minimuma napake glede na število členov asimptotske vrste.

Ustreznost referenčne funkcije¶

Za primerjavo naših vrst z referenčno vrednostjo lahko uporabimo funkciji mpmath.airyai() in mpmath.airybi() iz Pythonskega paketa mpmath. Ker bo to naša osnova za nadaljevanje, je smiselno preveriti, če ti dve funkciji dejansko delujeta pravilno oz. dovolj natančno glede na definicijo Airyjevih funkcij.

Z uporabo programa Wolfram Mathematica numerično integriramo definiciji funkcij $ \text{Ai} $ in $ \text{Bi} $ v nekaj testnih točkah s sledečimi nastavitvami:


xValues = Table[x, {x, -5, 15, 1}]; 
results = {};

integrandA[t_, x_] := Cos[t^3/3 + x*t];
integrandB[t_, x_] := Exp[-t^3/3 + x*t] + Sin[t^3/3 + x*t];

For[i = 1, i <= Length[xValues], i++, x = xValues[[i]];
 resultA = 
  NIntegrate[integrandA[t, x], {t, 0, Infinity}, PrecisionGoal -> 22, 
   AccuracyGoal -> 22, WorkingPrecision -> 60];
 resultB = 
  NIntegrate[integrandB[t, x], {t, 0, Infinity}, PrecisionGoal -> 22, 
   AccuracyGoal -> 22, WorkingPrecision -> 60];
 results = AppendTo[results, {x, resultA / Pi, resultB / Pi}];]

xValues = Table[x, {x, 15, 22, 1}]; 

For[i = 1, i <= Length[xValues], i++, x = xValues[[i]];
 resultA = 
  NIntegrate[integrandA[t, x], {t, 0, Infinity}, PrecisionGoal -> 20, 
   AccuracyGoal -> 20, WorkingPrecision -> 60];
 resultB = 
  NIntegrate[integrandB[t, x], {t, 0, Infinity}, PrecisionGoal -> 30, 
   AccuracyGoal -> 30, WorkingPrecision -> 60];
 results = AppendTo[results, {x, resultA / Pi, resultB / Pi}];]

Rezultate primarjamo s funkcijama mpmath.airyai() in mpmath.airybi(), ki ju prav tako nastavimo na 60 decimalk.

Odstopanje mpmath od numeričnega integrala: Primerjanih skupno 27 točk. Pri $x=15$ povečamo Precision in AccuracyGoal. Izkaže se, da je v Mathematici dobro targetirati PrecisionGoal in AccuracyGoal vsaj 5 redov višje, saj napaka v primerjavi pada.

Iz grafa razberemo, da mpmath.airyai() in mpmath.airybi() odstopada od dovolj natančno izračunanega integrala za manj kot $10^{-10}$. Torej lahko za primerjavo napak uporabimo kar ti dve funkciji, saj delujeta bistveno hitreje od numerične integracije.

Maclaurinov približek¶

Za majhne $x$ lahko funkciji $\text{Ai}$ in $\text{Bi}$ izrazimo z Maclaurinovima vrstama

$$ \text{Ai}(x) = \alpha f(x) - \beta g(x) $$

$$ \text{Bi}(x) = \sqrt{3}\, \Bigl[\alpha f (x) + \beta g(x) \Bigr] $$ kjer v $x=0$ velja $\alpha = \text{Ai}(0) \approx 0.3550...$ in $\beta = -\text{Ai}'(0) \approx 0.2588...$ Vrsti za $f$ in $g$ zapišemo rekurzivno in $n$-ti člen izračunamo kot produkt prejšnjega člena in količnika:

$$ f_{n} = f_{n-1} \cdot \frac{f_n}{f_{n-1}} = f_{n-1} \cdot \frac{x^3}{(3n - 1)\, 3n}, \quad f_0 = 1$$ $$ g_{n} = g_{n-1} \cdot \frac{g_n}{g_{n-1}} = g_{n-1} \cdot \frac{x^3}{(3n + 1)\, 3n}, \quad g_0 = x$$

Absolutna napaka: Iz grafov lahko opazimo, da absolutna napaka narašča z večanjem absolutne vrednosti $x$. Tik preden bi presegla mejo $10^{-10}$, se vrsti doda še en člen, ki napako zopet zmanjša.

Relativna napaka: Tudi tukaj opazimo podobno stvar. Tik preden bi relativna napaka presegla mejo $10^{-10}$, se vrsti doda še en člen, ki napako zopet zmanjša.

Asimptotski približek¶

Asimptotska vrsta lahko divergira ali konvergira. Napaka, ki jo naredimo pri krajšanju vrste za fiksen $x$ pri redu $n$ se ne zmanjšuje z večanjem reda $n$, kar pomeni, da se točki $x < \infty$ ne moremo poljubno približati z večanjem števila členov. Napaka se zmanjšuje proti nič šele, ko pri fiksnem $n$ povečujemo $x \rightarrow \infty$.

Za velike vrednosti $|x|$ Airyjevi funkciji aproksimiramo z njunima asimptotskima razvojema. Z novo spremenljivko $\xi=\frac{2}{3} |x|^{3/2}$ in asimptotskimi vrstami

$$ L(z) \sim \sum_{s=0}^\infty \frac{u_s}{z^s},\quad P(z) \sim \sum_{s=0}^\infty (-1)^s \frac{u_{2s}}{z^{2 s}},\quad Q(z) \sim \sum_{s=0}^\infty (-1)^s \frac{u_{2s+1}}{z^{2 s+1}}, $$

s koeficienti

$$ \begin{equation*} u_s = \frac{ \Gamma(3s + \frac{1}{2})} {54^s s!\, \Gamma(s + \frac{1}{2}) }. \end{equation*} $$

Za velike pozitivne $x$ izrazimo $$ \begin{equation*} \text{Ai}(x)\sim \frac{\mathrm{e}^{-\xi}}{2\sqrt{\pi} x^{1/4}} \, L(-\xi) \>, \quad \text{Bi}(x)\sim \frac{\mathrm{e}^{\xi}} { \sqrt{\pi} x^{1/4}} \, L(\xi)\>, \end{equation*} $$ za po absolutni vrednosti velike negativne $x$ pa $$ \begin{align*} \text{Ai}(x)&\sim \frac{1}{\sqrt{\pi} (-x)^{1/4}} \Bigl[ \phantom{-}\sin(\xi-\pi/4) \, Q(\xi) + \cos(\xi-\pi/4) \, P(\xi)\Bigr] \>, \\ \text{Bi}(x)&\sim \frac{1}{\sqrt{\pi} (-x)^{1/4}} \Bigl[ - \sin(\xi-\pi/4) \, P(\xi) + \cos(\xi-\pi/4) \, Q(\xi)\Bigr]\>. \end{align*} $$

Funkcije L, P in Q ponovno implementiramo rekurzivno.

Napaka za različno število členov: Asimptotskemu razvoju v našem primeru pada absolutna napaka samo do nekega končnega števila členov. Zato v nadaljevanju naloge jemljemo nove člene le, če so le-ti po absolutni vrednosti manjši od naslednjega.

Za velike pozitivne $x$¶

Za velike pozitivne $x$ implementiramo funkcijo $L(z)$ rekurzivno kot vsoto členov $l_n $ za $ n \in \mathbb{N}$

$$ l_{n} = l_{n-1} \cdot \frac{l_n}{l_{n-1}} = l_{n-1} \cdot \frac{(3n + \frac{5}{3})(3n + \frac{1}{2})}{(n + 1)\, 18 \, z}, \quad l_0 = 1.$$

Izračunane vrednosti primerjamo z referenčno funkcijo. Pri računanju z rekurzivno formulo povečujemo število členov dokler ne dosežemo najmanjšega razdalje do reference.

Absolutna napaka: Podobno kot pri prejšnjem razvoju opazimo stopničastost pri padcih napake. Število členov za doseganje ustrezne absolutne napake pri $\text{Bi}$ narašča. Vendar so tudi vrednosti funkcije ekstremne; $\text{Bi}$ izjemno hitro narašča, $\text{Ai}$ pa postaja zelo majhna.

Relativna napaka: Zanimivo, za doseganje ustrezne relativne napake potrebujemo vse manj členov, saj pri večjih vrednostih $x$ manjšo vlogo igrajo decimalke. Ta argument je smiseln za $\text{Bi}$, funkcija $\text{Ai}$ pa ima že sama po sebi pri tako velikih $x$ vrednost blizu $10^{-9}$, kar je skoraj znotraj iskanega intervala.

Za po absolutni vrednosti velike negativne $x$¶

Za po absolutni vrednosti velike negativne $x$ implementiramo funkciji $P(z)$ in $Q(z)$ rekurzivno kot vsoto členov $p_n $ oz. $q_n$ za $ n \in \mathbb{N}$

$$ p_{n} = p_{n-1} \cdot \frac{p_n}{p_{n-1}} = p_{n-1} \cdot - \frac{(6n - \frac{1}{2})(6n - \frac{5}{2})(6n - \frac{7}{2})(6n - \frac{11}{2})}{(2n - 1)(2n)\, 18^2 \, z^2}, \quad p_0 = 1$$

$$ q_{n} = q_{n-1} \cdot \frac{q_n}{q_{n-1}} = q_{n-1} \cdot - \frac{(6n - \frac{1}{2})(6n - \frac{5}{2})(6n + \frac{1}{2})(6n + \frac{5}{2})}{(2n + 1)(2n)\, 18^2 \, z^2}, \quad q_0 = \frac{5}{72z}.$$

Izračunane vrednosti primerjamo z referenčno funkcijo. Pri računanju z rekurzivno formulo število členov povečujemo dokler ne dosežemo najmanjšega razdalje do reference.

Absolutna napaka: Absolutna napaka je ustrezna, število členov z manjšanjem $x$ pada.

Relativna napaka: Relativna napaka je ustrezna, število členov in tudi relativna napaka z manjšanjem $x$ pada.

Zlepek približkov¶

V nadaljevanju bomo naše vrste zlepili in tako dobili Python funkcije, ki lahko učinkovito izračunajo $\text{Ai}$ in $\text{Bi}$ v okviru absolutne in relativne napake za celotno realno os.

Število členov v zlepkih: Na grafih prikažemo potrebno število členov za doseganje absolutne in relativne napake manjše od $10^{-10}$ za zlepke $\text{Ai}$ in $\text{Bi}$. Vsi dobljeni zlepki so relativno učinkoviti (potrebujejo majhno število členov) z izjemo zgornjega desnega grafa. Tam imamo v bližini $x=4$ težave z doseganjem ciljane napake, saj število členov tako za Maclaurinovo in asimptotsko vrsto podivja (>500).

Dodatno: Ničle funkcij¶

Ničle funkcije $\text{Ai}$ pogosto srečamo v matematični analizi pri določitvi intervalov ničel specialnih funkcij in ortogonalnih polinomov ter v fiziki pri računu energijskih spektrov kvantnomehanskih sistemov.

V nadaljevanju poiščemo prvih sto ničel $\{a_s\}_{s=1}^{100}$ Airyjeve funkcije $\text{Ai}$ in prvih sto ničel $\{b_s\}_{s=1}^{100}$ funkcije $\text{Bi}$ pri $x<0$. Za izračun uporabimo knjižnico mpmath in vgrajeni funkciji airyaizero in airybizero. Dobljene vrednosti primerjamo s formulama

$$ \begin{equation*} a_s = - f \left( \frac{3\pi(4s-1)}{8} \right) \>, \qquad b_s = - f \left( \frac{3\pi(4s-3)}{8} \right) \>, \qquad s = 1,2,\ldots \>, \end{equation*} $$ kjer ima funkcija $f$ asimptotski razvoj $$ \begin{equation*} f(z) \sim z^{2/3} \left( 1 + \frac{5}{48} \, z^{-2} -\frac{5}{36} \, z^{-4} +\frac{77125}{82944} \, z^{-6} -\frac{108056875}{6967296} \, z^{-8} + \ldots\right) \>. \end{equation*} $$

Ničle Airyjevih funkcij: Iz rezultatov potrdimo, da zgornja asimptotska vrsta znotraj intervalov napak ustreza ničlam Airyjevih funkcij. Absolutna napaka se z večanjem $n$ manjša.

Zaključek¶

Tekom tega poročila nam je uspelo poiskati učinkovit način za izračun Airyjevih funkcij z absolutno napako manjšo od $10^{-10}$. Tipično število potrebovanih členov vrste je manjše od 25. Potencialne težave imamo le pri funkciji $\text{Bi}$ za velike $x$, saj tam število potrebnih členov izgleda kot da počasi narašča. Vendar v tem območju funkcija presega vrednosti $10^{100}$, tako da se nam doseganje absolutne napake manjše od $10^{-10}$ ne zdi pretirano praktično uporabno. Podoben uspeh smo dosegli tudi za relativno napako manjšo od $10^{-10}$, kjer sta funkciji učinkoviti na celotni realni osi z izjemo okolice $x=4$ pri funkciji $\text{Ai}$. Tam število potrebovanih členov strmo naraste (na več kot 500) za asimptotsko in Maclaurinovo vrsto, posledično v tej okolici za majhno število členov ciljane relativne napake ne dosežemo.

Za konec smo izračunali še prvih 100 ničel obeh Airyjevih funkcij in jih primerjali z asimptotskim razvojem. Pričakovano, se absolutna napaka manjša z manjšanjem $x$ (negativna vrednost) oz. večanjem $n$ (tj. indeksa, ki označuje ničle na negativni polosi).

Luka Skeledžija, Github source 🔗, 2023

01 - Airyjevi funkciji¶